library(tidyverse)
library(sf)
library(terra)
library(ggpubr)
ggplot2::theme_set(theme_classic())
This document visualizes the data of vegetation cover by functional type that we use in subsequent analysis at each step of the data-processing workflow. These maps serve both as a reference and a check that the workflow is performing as intended.
For the sake of simplified comparison, I show data from 2016, 2020, and 2022, which are years for which we have easy to access rasters of FIA data to compare to.
## Reading layer `vegCompPoints' from data source
## `/Users/astears/Documents/Dropbox_static/Work/NAU_USGS_postdoc/cleanPED/PED_vegClimModels/Data_processed/CoverData/DataForAnalysisPoints'
## using driver `ESRI Shapefile'
## Simple feature collection with 1266274 features and 30 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -124.7004 ymin: 24.68798 xmax: -67.10236 ymax: 49.3489
## Geodetic CRS: WGS 84
# Read in data
dat_2 <- readRDS(file = "../../../Data_processed/CoverData/dataForAnalysis_fireRemoved.rds")
# trim to be only later than 2000, which is what we want to use for further analysis
dat_2 <- dat_2 %>%
filter(Year >= 2000)
## make figure for RAP
# get data just for RAP
RAP_2 <- dat_2 %>%
filter(Source == "RAP") %>%
pivot_longer(cols = c(ShrbCvr, TtlTrCv, TtlHrbC, BrGrndC#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
),
names_to = "coverType",
values_to = "coverValues"
) %>%
filter(!is.na(coverValues))
ggarrange(ggplot(RAP_2 %>%
filter(Year %in% c(2016, 2020, 2022))) +
facet_grid(rows = vars(coverType), cols = vars(Year)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("Filtered Burned Areas - RAP"),
ggplot(RAP_2) +
facet_grid(rows = vars(coverType)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("Filtered Burned Areas, \nall years - RAP "),
ncol = 2, align = "h",
common.legend = TRUE,
widths = c(.7, .3)
)
# Read in data
dat_3 <- readRDS(file = "../../../Data_processed/CoverData/data_beforeSpatialAveraging_sampledLANDFIRE.rds")
# trim to be only later than 2000, which is what we want to use for further analysis
dat_3 <- dat_3 %>%
filter(Year >= 2000)
## make figure for RAP
# get data just for RAP
RAP_3 <- dat_3 %>%
filter(Source == "RAP") %>%
pivot_longer(cols = c(#ShrbCvr, TtlTrCv, TtlHrbC, BrGrndC
ShrubCover, TotalTreeCover, TotalHerbaceousCover, BareGroundCover#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
),
names_to = "coverType",
values_to = "coverValues"
) %>%
filter(!is.na(coverValues))
ggarrange(ggplot(RAP_3 %>%
filter(Year %in% c(2016, 2020, 2022))) +
facet_grid(rows = vars(coverType), cols = vars(Year)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("After 'coarsifying' LANDFIRE data - RAP"),
ggplot(RAP_3) +
facet_grid(rows = vars(coverType)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("After 'coarsifying' LANDFIRE data, \nall years - RAP "),
ncol = 2, align = "h",
common.legend = TRUE,
widths = c(.7, .3)
)
## make figure for LDC
# get data just for LDC
LDC_1 <- dat_1 %>%
st_zm() %>%
filter(Source == "LDC") %>%
#st_drop_geometry() %>%
pivot_longer(cols = c(ShrbCvr, TtlTrCv, TtlHrbC, BrGrndC#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
),
names_to = "coverType",
values_to = "coverValues"
) %>%
filter(!is.na(coverValues))
#filter(coverType == "ShrbCvr" & Year == 2016) %>%
ggarrange(ggplot(LDC_1 %>%
filter(Year %in% c(2016, 2020, 2022))) +
facet_grid(rows = vars(coverType), cols = vars(Year)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("Initial Data - AIM"),
ggplot(LDC_1) +
facet_grid(rows = vars(coverType)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("Initial Data, \nall years - AIM "),
ncol = 2, align = "h",
common.legend = TRUE,
widths = c(.7, .3)
)
## make figure for LDC
# get data just for LDC
AIM_2 <- dat_2 %>%
filter(Source == "LDC") %>%
pivot_longer(cols = c(ShrbCvr, TtlTrCv, TtlHrbC, BrGrndC#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
),
names_to = "coverType",
values_to = "coverValues"
) %>%
filter(!is.na(coverValues))
ggarrange(ggplot(AIM_2 %>%
filter(Year %in% c(2016, 2020, 2022))) +
facet_grid(rows = vars(coverType), cols = vars(Year)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("Filtered Burned Areas - AIM"),
ggplot(AIM_2) +
facet_grid(rows = vars(coverType)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("Filtered Burned Areas, \nall years - AIM "),
ncol = 2, align = "h",
common.legend = TRUE,
widths = c(.7, .3)
)
## make figure for AIM
# get data just for AIM
AIM_3 <- dat_3 %>%
filter(Source == "LDC") %>%
pivot_longer(cols = c(ShrubCover, TotalTreeCover, TotalHerbaceousCover, BareGroundCover#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
),
names_to = "coverType",
values_to = "coverValues"
) %>%
filter(!is.na(coverValues))
ggarrange(ggplot(AIM_3 %>%
filter(Year %in% c(2016, 2020, 2022))) +
facet_grid(rows = vars(coverType), cols = vars(Year)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("After 'coarsifying' LANDFIRE data - AIM"),
ggplot(AIM_3) +
facet_grid(rows = vars(coverType)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("After 'coarsifying' LANDFIRE data, \nall years - AIM "),
ncol = 2, align = "h",
common.legend = TRUE,
widths = c(.7, .3)
)
## make figure for FIA
# get data just for FIA
FIA_1 <- dat_1 %>%
st_zm() %>%
filter(Source == "FIA") %>%
#st_drop_geometry() %>%
pivot_longer(cols = c(ShrbCvr, TtlTrCv, TtlHrbC, BrGrndC#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
),
names_to = "coverType",
values_to = "coverValues"
) %>%
filter(!is.na(coverValues))
ggarrange(ggplot(FIA_1 %>%
filter(Year %in% c(2016, 2020, 2022))) +
facet_grid(rows = vars(coverType), cols = vars(Year)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("Initial Data - FIA"),
ggplot(FIA_1) +
facet_grid(rows = vars(coverType)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("Initial Data, \nall years - FIA "),
ncol = 2, align = "h",
common.legend = TRUE,
widths = c(.7, .3)
)
## make figure for FIA
# get data just for FIA
FIA_2 <- dat_2 %>%
filter(Source == "FIA") %>%
pivot_longer(cols = c(ShrbCvr, TtlTrCv, TtlHrbC, BrGrndC#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
),
names_to = "coverType",
values_to = "coverValues"
) %>%
filter(!is.na(coverValues))
ggarrange(ggplot(FIA_2 %>%
filter(Year %in% c(2016, 2020, 2022))) +
facet_grid(rows = vars(coverType), cols = vars(Year)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("Filtered Burned Areas - FIA"),
ggplot(FIA_2) +
facet_grid(rows = vars(coverType)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("Filtered Burned Areas, \nall years - FIA "),
ncol = 2, align = "h",
common.legend = TRUE,
widths = c(.7, .3)
)
## make figure for FIA
# get data just for FIA
FIA_3 <- dat_3 %>%
filter(Source == "FIA") %>%
pivot_longer(cols = c(ShrubCover, TotalTreeCover, TotalHerbaceousCover, BareGroundCover#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
),
names_to = "coverType",
values_to = "coverValues"
) %>%
filter(!is.na(coverValues))
ggarrange(ggplot(FIA_3 %>%
filter(Year %in% c(2016, 2020, 2022))) +
facet_grid(rows = vars(coverType), cols = vars(Year)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("After 'coarsifying' LANDFIRE data - FIA"),
ggplot(FIA_3) +
facet_grid(rows = vars(coverType)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("After 'coarsifying' LANDFIRE data, \nall years - FIA "),
ncol = 2, align = "h",
common.legend = TRUE,
widths = c(.7, .3)
)
## make figure for LANDFIRE
# get data just for LANDFIRE
LANDFIRE_1 <- dat_1 %>%
st_zm() %>%
filter(Source == "LANDFIRE") %>%
#st_drop_geometry() %>%
pivot_longer(cols = c(ShrbCvr, TtlTrCv, TtlHrbC, BrGrndC#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
),
names_to = "coverType",
values_to = "coverValues"
) %>%
filter(!is.na(coverValues))
ggarrange(ggplot(LANDFIRE_1 %>%
filter(Year %in% c(2003, 2007, 2015))) +
facet_grid(rows = vars(coverType), cols = vars(Year)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("Initial Data - LANDFIRE"),
ggplot(LANDFIRE_1) +
facet_grid(rows = vars(coverType)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("Initial Data, \nall years - LANDFIRE "),
ncol = 2, align = "h",
common.legend = TRUE,
widths = c(.7, .3)
)
## make figure for LANDFIRE
# get data just for LANDFIRE
LANDFIRE_2 <- dat_2 %>%
filter(Source == "LANDFIRE") %>%
pivot_longer(cols = c(ShrbCvr, TtlTrCv, TtlHrbC, BrGrndC#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
),
names_to = "coverType",
values_to = "coverValues"
) %>%
filter(!is.na(coverValues))
ggarrange(ggplot(LANDFIRE_2 %>%
filter(Year %in% c(2003, 2007, 2015))) +
facet_grid(rows = vars(coverType), cols = vars(Year)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("Filtered Burned Areas - LANDFIRE"),
ggplot(LANDFIRE_2) +
facet_grid(rows = vars(coverType)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("Filtered Burned Areas, \nall years - LANDFIRE "),
ncol = 2, align = "h",
common.legend = TRUE,
widths = c(.7, .3)
)
## make figure for LANDFIRE
# get data just for LANDFIRE
LANDFIRE_3 <- dat_3 %>%
filter(Source == "LANDFIRE") %>%
pivot_longer(cols = c(ShrubCover, TotalTreeCover, TotalHerbaceousCover, BareGroundCover#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
),
names_to = "coverType",
values_to = "coverValues"
) %>%
filter(!is.na(coverValues))
ggarrange(ggplot(LANDFIRE_3 %>%
filter(Year %in% c(2003, 2007, 2015))) +
facet_grid(rows = vars(coverType), cols = vars(Year)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("After 'coarsifying' LANDFIRE data - LANDFIRE"),
ggplot(LANDFIRE_3) +
facet_grid(rows = vars(coverType)) +
stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) +
scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"),
limits = c(0,100)) +
#geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) +
#geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
xlim(c(-130,-65)) +
ylim(c(25, 50)) +
ggtitle("After 'coarsifying' LANDFIRE data, \nall years - LANDFIRE "),
ncol = 2, align = "h",
common.legend = TRUE,
widths = c(.7, .3)
)
## After adding climate/weather data
## After adding ecoregion data